/*	Interactive Gravitational Simulator
	Project: IGS Command Line
	File: main.cpp
	Author: Mike Bantegui <mbante2@gmail.com>, Hofstra University 
	Copyright (C) 2012 - 2013 Mike Bantegui

	This file is part of the IGS software package for simulating
	the N-Body problem in real-time. This was developed as part
	of a departmental honors thesis project at Hofstra University
	in the Spring 2012 semester.

	IGS is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2 of the License, or
	(at your option) any later version.

	IGS is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
*/

#include "Body.h"
#include "Integrator.h"
#include "Model.h"
#include "NBodySystem.h"
#include "Octree.h"
#include "Simulator.h"
#include <iostream>

int main()
{
	using namespace IGS;

	int n = 2;
	Body *bodies = new Body[n];

	// Test Orbit: Circular orbit
	bodies[0].Mass = 1;
	bodies[0].Position = Vector(+1, 0, 0);
	bodies[0].Velocity = Vector(0, +0.5, 0);
	bodies[0].Acceleration.Zero();
	bodies[1].Mass = 1;
	bodies[1].Position = Vector(-1, 0, 0);
	bodies[1].Velocity = Vector(0, -0.5, 0);
	bodies[1].Acceleration.Zero();

	Simulator *sim = new Simulator(new NBodySystem(bodies, n), new EulerTrapezoid());
	std::shared_ptr<NBodySystem> system = sim->GetSystem();

	double pi = 2 * asin(1.0);
	double R = Distance(bodies[0].Position, bodies[1].Position);
	double T = 2 * pi * sqrt(R * R * R / (bodies[0].Mass + bodies[1].Mass));

	int steps = 1000;
	double dt = T / steps;
	double E0 = system->TotalEnergy();

	// Print position of all bodies
	std::cout << "Positions:" << std::endl;
	for (int i = 0; i < n; i++)
		std::cout << "Body " << i << ": (" << bodies[i].Position << ")" << std::endl;
	std::cout << "N = " << system->Size() << std::endl;
	std::cout << "KE = " << system->KineticEnergy() << std::endl;
	std::cout << "PE = " << system->PotentialEnergy() << std::endl;
	std::cout << "E = " << E0 << std::endl;
	std::cout << "Center of Mass = (" << system->CenterOfMass() << ")" << std::endl;
	std::cout << "Center of Momentum = (" << system->CenterOfMomentum() << ")" << std::endl;
	std::cout << "-----------------------------------" << std::endl;

	std::cout << "Hit enter to continue" << std::endl;
	std::cin.get();

	for (int step = 0; step < steps; step++)
	{
		sim->Evolve(dt);
		double E = system->TotalEnergy();
		std::cout << "E = " << E;
		std::cout << ", % Error = " << 1 - E / E0 << std::endl;
		std::cout << "Center of Mass = (" << system->CenterOfMass() << ")" << std::endl;
		std::cout << "Center of Momentum = (" << system->CenterOfMomentum() << ")" << std::endl;
	}

	for (int i = 0; i < n; i++)
		std::cout << "(" << bodies[i].Position << ")" << std::endl;
	std::cout << "E = " << system->TotalEnergy() << std::endl << std::endl;

	NBodySystem::Write(*system, "Model.txt");

	// Wait for user to step
	std::cin.get();

	return 0;
}